#include <iostream>
#include <string>
#include <fstream>

#include "VortSystem.h"
#include "../tinyxml/tinyxml.h"

using namespace VortPan;
using namespace std;


void WritePoints(std::vector<Vector3D>& Stream, std::string FName){
	ofstream of(FName.c_str());
	
	for (std::vector<Vector3D>::iterator it=Stream.begin(); it!=Stream.end(); it++) {
		of<<it->v[0]<<" "<<it->v[1]<<" "<<it->v[2]<<endl;
	}
	
}

Vector3D ReadVector(TiXmlElement *pEl){
	Vector3D v=V0();
	
	v.v[0]=atof(pEl->FirstChildElement("x")->GetText());
	v.v[1]=atof(pEl->FirstChildElement("y")->GetText());
	v.v[2]=atof(pEl->FirstChildElement("z")->GetText());

	if(pEl->FirstChildElement("Normalize")){
		v=Normalize(v);
	}
	
	
	return v;
}



int main (int argc, char * const argv[]) {
	string Name;
	
	if(argc<2){
		Name="case.xml";
	} else {
		Name=argv[1];
	}

	TiXmlDocument Doc;
	assert(Doc.LoadFile(Name.c_str()));

	TiXmlElement *pRoot;
	assert(pRoot=Doc.FirstChildElement());
	
	VortMesh Mesh;
	
	double WeldingTolerance=1e-5;
	
	TiXmlElement *pElem;
	if(pElem=pRoot->FirstChildElement("WeldingTolerance")){
		WeldingTolerance=atof(pElem->GetText());
	}
	
	re Scale=1e-0;
	
	TiXmlElement *pMeshes;
	
	assert(pMeshes=pRoot->FirstChildElement("Meshes"));
	
	TiXmlElement *pEl=NULL;
	string MeshFile;
	bool IsWake;
	
	
	for(pEl=pMeshes->FirstChildElement("Mesh"); pEl!=NULL; pEl=pEl->NextSiblingElement()) {
		VortMesh LoadMesh;

		MeshFile=pEl->GetText();
		Scale=1.0;

		TiXmlElement *pParamElement;

		pParamElement=pEl->FirstChildElement("Scale");
		if (pParamElement!=NULL) {
			Scale=atof(pParamElement->GetText());
		}

		IsWake=false;
		if(pEl->FirstChildElement("Wake")){
			IsWake=true;
		}
		
		LoadMesh.LoadPLY(MeshFile, IsWake, Scale);
				
		Mesh.AddMesh(LoadMesh);
	}
	
	
	
	Mesh.Weld(WeldingTolerance);

	

	

	
	
	Vector3D Velocity=Normalize(e1());

	assert(pElem=pRoot->FirstChildElement("Velocity"));
	
	Velocity=ReadVector(pElem);
	
	
	int NIter=1;
	re RelaxationFactor=1;
	re VortexCoreSize=0.1;

	if(pElem=pRoot->FirstChildElement("WakeRelaxation")){
		
		Vector3D WakeDir=Normalize(Velocity);

		TiXmlElement *pParamElement;
		
		if(pParamElement=pElem->FirstChildElement("ReplaceWake")){
			
			if (pEl=pParamElement->FirstChildElement("InitialWakeDirection")){
				WakeDir=ReadVector(pEl);
			}
			assert(pEl=pParamElement->FirstChildElement("StartSegmentLength"));
			re StartSegmentLength=atof(pEl->GetText());
			
			assert(pEl=pParamElement->FirstChildElement("EndDistance"));
			re EndDistance=atof(pEl->GetText());
			
			assert(pEl=pParamElement->FirstChildElement("NextElementRatio"));
			re NextElementRatio=atof(pEl->GetText());
			
			Mesh.ReplaceAllWakePolys(WakeDir, StartSegmentLength, EndDistance, NextElementRatio, WeldingTolerance);
		}
		
		assert(pEl=pElem->FirstChildElement("NumberOfIterations"));
		NIter=atof(pEl->GetText());

		assert(pEl=pElem->FirstChildElement("RelaxationFactor"));
		RelaxationFactor=atof(pEl->GetText());

		assert(pEl=pElem->FirstChildElement("VortexCoreSize"));
		VortexCoreSize=atof(pEl->GetText());
		
	}
	
	VortSystem System(Mesh);
	
	
	TiXmlElement *pPlanes;

	if (pPlanes=pRoot->FirstChildElement("SymmetryPlanes")) {
		for(pEl=pPlanes->FirstChildElement("Plane"); pEl!=NULL; pEl=pEl->NextSiblingElement()) {
			Vector3D Normal,Position;
			Position=V0();
			Normal=e3();
			TiXmlElement *pParamElement;
			if(pParamElement=pEl->FirstChildElement("Position")){
				Position=ReadVector(pParamElement);
			}
			if(pParamElement=pEl->FirstChildElement("Normal")){
				Normal=ReadVector(pParamElement);
			}
			VortSymmetryPlane Plane(Position,Normal);
			System.AddSymmetry(Plane);
		}
	}
	
	System.CheckMesh(WeldingTolerance);
	
	
	System.SetVelocity(Velocity);

	

	for (int i=0; i<NIter; i++) {
		System.SolveSystem();
		
		System.RelaxWake(RelaxationFactor,VortexCoreSize);		
	}

//	System.GetMesh().ExportMathematicaList("MathWing.dat");
//	System.GetMesh().ExportCenterPoints("MathCenter.dat");
	
/*	std::vector<Vector3D> Stream;
	Vector3D R0=0.3*e2()-0.8*e1()-0.05*e3();
	System.StreamLine(R0, 0.01, 10, Stream);	
	WritePoints(Stream, "TestPoints.dat");
	System.StreamLineRegular(R0, 0.01, 10, Stream, 0.01);	
	WritePoints(Stream, "TestPointsRegular.dat");*/
	
	Memory Pressures,Gradients;
	
	Vector3D F,M;
	
	System.CalculateDerivativePressures(Pressures,&F,&M);
	
	cout<<"Force= "<<F.v[0]<<", "<<F.v[1]<<", "<<F.v[2]<<",  Moment= "<<M.v[0]<<", "<<M.v[1]<<", "<<M.v[2]<<endl;
	
	string PLYName="finalmesh.ply";
	
	if (pElem=pRoot->FirstChildElement("FinalMeshName")) {
		PLYName=pElem->GetText();
	}
	
	System.GetMesh().SavePLY(PLYName);
	
	System.CalculateAbsolutePressureGradient(Pressures,Gradients);

	if (pElem=pRoot->FirstChildElement("PressuresName")) {
		string PName=pElem->GetText();
		System.ExportData(PName,Pressures);
	}

	if (pElem=pRoot->FirstChildElement("PressureGradientsName")) {
		string PName=pElem->GetText();
		System.ExportData(PName,Gradients);
	}
	
	//System.ExportData("TestPressures.dat",Pressures);
	//System.ExportData("TestGradients.dat",Gradients);
	//System.ExportIsWake("TestWake.dat");
	//System.ExportClC("TestClC.dat");
	
//	cout<<"Saving="<<Mesh.ExportMathematicaList("TestWing.dat")<<endl;
//	cout<<"BadEdgesExport="<<Mesh.ExportBadEdges("BadEdges.dat")<<endl;
	
	
}
